以复现图表的方式来学习一篇文章
分享是一种态度
1.文章简介
本次学习的2019年7月发表于Circulation.文章题目为：Single-Cell Analysis of the Normal Mouse Aorta Reveals Functionally Distinct Endothelial Cell Populations
”
2. 文章思路很直接套路
做Normal Aorta的 Single-Cell,对其进行分群，其中有一个亚群叫EC 重点关注亚群EC，对EC群再分群，做gene set enrichment 分析，做marker gene - 表达的小提琴图，heatmap图 实验验证EC群的特性。比如不同的diet之后(心血管疾病肯定和饮食之类的因素有关)，做不同condition 下的single-Cell，看EC群基因变化，然后就和疾病关系靠一靠
作业复现的图表是Figure1 的1B 和1C
3. 文件下载和解读
主要用Seurat V3. Broad Institute的single cell portal上面存放了他们发表的single-cell文章的数据 这里一共是3个文件
Meta_DATA_Chow_12PCs_outfile.gz Cluster_File_Chow_12PCs_outfile.gz Seurat_Chow_12PCs_outfile.mtx 通过命名可以看出，文件均为跑过了PCA，tSNE分群后的输出数据，也就是说，这次的任务是非常下游的可视化过程。
读取文件1和2
安装加载包，我用的是Seurat v3
library(data.table)
library(R.utils)
library(Seurat)
library(ggplot2)
读取文件1和文件2
cluster_file <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/Cluster_File_Chow_12PCs_outfile.gz',
header = T,sep = '\t',data.table = F)
meta_data <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/META_DATA_Chow_12PCs_outfile.gz',
header = T, sep = '\t',data.table = F)
这里看一下两个文件。cluster_file显示的是tSNE_x和tSNE_y的坐标。这也是我们画fig1B tSNE图需要的。这个nametype列名是样品名 meta_data显示了分群信息。在Jimmy大神前面单细胞的讲座中，了解到需要将这些数据的行名(rownames)改为样品名。
head(cluster_file[1:3,1:3])
## NAME X Y
## 1 TYPE numeric numeric
## 2 AAACCTGAGACGCTTT-1 -19.0225914077646 0.924488198164653
## 3 AAACCTGCACTGTGTA-1 -26.7488891822262 -0.784745741413262
head(meta_data[1:3,1:3])
## NAME Cluster Sub-Cluster
## 1 TYPE group group
## 2 AAACCTGAGACGCTTT-1 VSMC VSMC
## 3 AAACCTGCACTGTGTA-1 VSMC VSMC
读取文件3 我在读取.mtx的文件上，入了好大的坑。看到后缀是.mtx文件，我的第一个反应是用readMM()读。但是，持续报错，error msg说这个文件不是稀疏矩阵。后来还是和前面一样用fread()读的，居然读出来了。然后我看了一下，是每个基因在每个样品中的表达，在数值上也有稀疏矩阵的.这个符号。读入后，我用了typeof()看这个.mxt，它显示的是list.这一部分我至今困惑。言归正传，为了后续的分析，需要将读入的.mtx文件的行名改成基因名。
mxt <- fread(file = '/Users/xiaoming/Downloads/2019_circulation/data/Seurat_Chow_12PCs_outfile.mtx',header = T, sep = '\t',data.table = F)
head(mxt[1:3,1:3])
## GENE AAACCTGAGACGCTTT-1 AAACCTGCACTGTGTA-1
## 1 0610007P14Rik 0.7264328 0.0000000
## 2 0610009B22Rik 0.0000000 0.6094115
## 3 0610009L18Rik 0.7264328 0.0000000
4. 重现fig1B
由于作者已经给出了tSNE的x，y轴坐标，这一步就特别简单。直接把有坐标信息的Cluster_File_Chow_12PCs_outfile.gz文件和有分群信息的Meta_DATA_Chow_12PCs_outfile.gz一整合，生成一个既有坐标信息，又有分群信息的文件（该文件叫dat),然后用ggplot画图就可以。
# 整合坐标信息和分群信息的文件
phe=data.frame(cell=meta_data$NAME[-1],
cluster =meta_data$Cluster[-1])
cluster_file_tSNE <- data.frame(cell = cluster_file$NAME[-1],
tSNE_1 =cluster_file$X[-1],
tSNE_2 = cluster_file$Y[-1],
avg_intensity =meta_data$`Average Intensity`[-1])
dat=merge(cluster_file_tSNE,phe,by.x = "cell",by.y="cell")
dat$tSNE_1 <- as.numeric(dat$tSNE_1)
dat$tSNE_2 <- as.numeric(dat$tSNE_2)
#画出图1B
dat_tsne =ggplot(dat,aes(x=tSNE_1,y=tSNE_2,color=cluster))+geom_point(size=0.95)
5.重现fig1C
fig1C的主体思路是找到每个群的marker gene,然后用Seurat v3自带Dotplot()可视化。marker gene的找法用Seurat自带的FindMarkers().
首先FindMarkers()里面的一个输入对象是object，这里需要使用CreateSeuratObject用.mtx这个文件创建一个对象。
# 创建Seurat对象
mxt_obj <- CreateSeuratObject(counts = mxt, project = "mxt", min.cells = 3, min.features = 200)
找marker gene需要分群信息。这里可以深刻体会AddMetaData()的意义。通过AddMetaData(),给CreateSeuratObject创建的对象提供更多的信息。作者已经在metadata文件里面提供了分群，那么通过AddMetaData()，就直接可以把分群信息赋予创建的对象，不用再通过单细胞分析全流程FindCluster()等函数计算。
另外，需要注意的是，如果想用FindMarkers()这个函数找差异基因，它认定的分群需要用数字表示，不能识别字符串比如’VSMC’这种分群，所以得把作者提供的meta-data里面的分群信息改成数字1，2表示
# 将分群信息改成数字形式
meta_data_3 <- meta_data
meta_data_3$Cluster <- as.numeric(meta_data_3$Cluster)
head(meta_data_3)
## Cluster
## AAACCTGAGACGCTTT-1 6
## AAACCTGCACTGTGTA-1 6
## AAACGGGCACGGATAG-1 6
## AAACGGGGTGTTCTTT-1 6
## AAAGATGAGGCTCATT-1 6
## AAAGATGCAAGACACG-1 6
# 创建的对象中使用AddMetaData()添加分群信息
mxt_obj <- AddMetaData(object = mxt_obj,
metadata = meta_data_3,
col.name = 'Cluster')
进入FindMarkers()这个函数。注意两个参数，group.by 和ident.1
ident.1参数指明你要找的marker gene 所在的群。前面说了，ident.1参数只认数字代表的分群，不认识字符串。所以，在meta-data里面，将分群信息改为数字
group.by是划分的标准，比如，你找不同分群的marker gene，分群信息就是在meta_data里面的哪一列保存，列名就是group.by的参数值。在这里，我们这个分群信息保存在列名为"Cluster"的列下面
# find marker gene for clsuter VSMC。VSMC对应的群的数字是6
markers_VSMC <- FindMarkers(mxt_obj,
group.by = 'Cluster',
ident.1 = 6,
test.use = "wilcox",
min.pct = 0.25)
文中提到，marker gene 的定义是 genes with the highest differential expression relative to all other genes. 所以得到的结果就按照avg_logFC排序，取前5。
注明：我找的marker gene和文章的做对比。就EC群的Ccl21a对不上。作者找的叫Ccl21a,我找了个，叫Cldn5.两不是一个基因
全部marker gene得到后，用seurat 自带的dotplot()画图
DotPlot(object = mxt_meta2, features = Marker_gene,group.by = 'Cluster')+RotatedAxis()
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
生信爆款入门-第7期（线上直播4周，马拉松式陪伴，带你入门）你的生物信息入门课
数据挖掘学习班第5期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课
生信技能树的2019年终总结 你的生物信息成长宝藏
看完记得顺手点个“在看”哦！
长按扫码可关注